!************************************************************
! this subrouitne calculate direct wall point of isentropic
! supersonic flow
!************************************************************
! Created by : B. G.
! Date       : 2015-08-12
! Revised    :
!************************************************************
! Nomenclature:
!   yp~rho0   temporary properties, p-c+, 0-streanline
!   lp~T02    coefficients of equations
!   *4f       previous loop properties of point 4
!   i         loop indicator   
!   
!************************************************************
! Warning :
!   1. angle in rad, not degree
!************************************************************

subroutine DWall()
   !          4()
   ! 3
   !      2      
   use VariableDef
   implicit none
   ! given variables
   ! real*8::x2,y2,V2,theta2,p2,rho2,x3,y3,V3,theta3,p3,rho3
   ! temporary variables
   real*8::yp,pp,rhop,Vp,alphap,thetap,ap,Mp
   real*8::p00,V0,rho0
   real*8::lp,Qp,Sp,Tp,R0,A0,T01,T02
   real*8::x4f,y4f,p4f,rho4f,V4f,theta4f
   integer::i
   ! calculate variables
   ! real*8::x4,y4,V4,theta4,p4,rho4,T4,M4
   
   yp=y2
   Vp=V2
   thetap=theta2
   pp=p2
   rhop=rho2
   
   p4=p3
   V4=V3
   rho4=rho3
   
   x4f=0.0
   y4f=0.0
   p4f=0.0
   rho4f=0.0
   V4f=0.0
   theta4f=0.0
   
   do i=1,icor
      ! determine coefficients for l-c-line
      call Thermo(pp,rhop,Vp,ap,T4,Mp)
      alphap=asin(1/Mp)
      lp=tan(thetap+alphap)
      Qp=sqrt(Mp**2-1)/(rhop*Vp**2)
      Sp=delta*sin(thetap)/(yp*Mp*cos(thetap+alphap))
      ! determine x4, y4, theta4
      call Wall24(lp)
      ! determine Tp
      Tp=-Sp*(x4-x2)+Qp*p2+theta2
      ! determine coefficients for streamline
      p00=(p3+p4)/2
      V0=(V3+V4)/2
      rho0=(rho3+rho4)/2
      call Thermo(p00,rho0,V0,ap,T4,M4)
      R0=rho0*V0
      A0=ap**2
      T01=R0*V3+p3
      T02=p3-A0*rho3
      ! determine p4 V4 rho4
      p4=(Tp-theta4)/Qp
      V4=(-p4+T01)/R0
      rho4=(p4-T02)/A0
      ! justify if point 4 arrived
      if ((abs(x4-x4f) <= el) .and. (abs(y4-y4f) <= el) .and.&
         &(abs(p4-p4f) <= ep*p4f) .and. (abs(rho4-rho4f) <= erho*rho4f) .and.&
         &(abs(V4-V4f) <= eV*V4f) .and. (abs(theta4-theta4f) <= etheta*theta4f)) exit
      
      ! reset c+ and c- values
      yp=(y2+y4)/2
      pp=(p2+p4)/2
      rhop=(rho2+rho4)/2
      Vp=(V2+V4)/2
      thetap=(theta2+theta4)/2
      ! store as front loop values
      x4f=x4
      y4f=y4
      p4f=p4
      rho4f=rho4
      V4f=V4
      theta4f=theta4
   end do
   call Thermo(p4,rho4,V4,ap,T4,M4)
   Vv4=V4*sin(theta4)
   Vu4=V4*cos(theta4)

end subroutine DWall
